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A fundamental problem in remote sensing and radiative transfer simulations involving ice 
clouds is the ability to compute accurate optical properties for individual ice particles. While 
relatively simple and intuitively appealing, the conventional geometric-optics method 
(CGOM) is used frequently for the solution of light scattering by ice crystals. Due to the 
approximations in the ray-tracing technique, the CGOM accuracy is not well quantified. The 
result is that the uncertainties are introduced that can impact many applications. Improve- 
ments in the Invariant Imbedding T-matrix method (II-TM) and the Improved Geometric- 
Optics Method (IGOM) provide a mechanism to assess the aforementioned uncertainties. 
The results computed by the II-TM + IGOM are considered as a benchmark because the II- 
TM solves Maxwell's equations from first principles and is applicable to particle size 
parameters ranging into the domain at which the IGOM has reasonable accuracy. To assess 
the uncertainties with the CGOM in remote sensing and radiative transfer simulations, two 
independent optical property datasets of hexagonal columns are developed for sensitivity 
studies by using the CGOM and the II-TM + IGOM, respectively. Ice cloud bulk optical 
properties obtained from the two datasets are compared and subsequently applied to 
retrieve the optical thickness and effective diameter from Moderate Resolution Imaging 
Spectroradiometer (MODIS) measurements. Additionally, the bulk optical properties are 
tested in broadband radiative transfer (RT) simulations using the general circulation model 
(GCM) version of the Rapid Radiative Transfer Model (RRTMG) that is adopted in the 
National Center for Atmospheric Research (NCAR) Community Atmosphere Model (CAM, 
version 5.1). For MODIS retrievals, the mean bias of uncertainties of applying the CGOM in 
shortwave bands (0.86 and 2.13 pm) can be up to 5% in the optical thickness and as high as 
20% in the effective diameter, depending on cloud optical thickness and effective diameter. 
In the MODIS infrared window bands centered at 8.5, 11, and 12 pm, biases in the optical 
thickness and effective diameter are up to 12% and 10%, respectively. The CGOM-based 
simulation errors in ice cloud radiative forcing calculations are on the order of 10 W m " 2 . 
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1. Introduction 

Ice clouds are known to play an important role in 
regulating atmospheric radiation through interaction with 
both solar and infrared (IR) radiation fields [1,2], Theore- 
tical modeling of ice cloud radiative effects is indispensa- 
ble in many atmospheric applications. The modeling 
approach typically begins with the single scattering of 
light by individual ice particles of different habits (shapes) 
and sizes contained in a scattering volume element, 
followed by the multiple-scattering calculations (i.e„ the 
radiative transfer, or RT) involving ice clouds. For example, 
inferring the microphysical properties of ice clouds from 
observations by ground- and satellite-based instruments 
requires comparing the instrument observations to the 
model-simulated results obtained from the single- 
scattering properties and relevant multiple-scattering 
simulations [3-7]. To understand the ice cloud impact in 
the climate models, an accurate representation of ice cloud 
bulk-scattering properties is a critical component required 
by RT models to quantify the cloud radiative effect (CRE), 
such as radiative forcing of natural cirrus clouds [8,9] and 
anthropogenic contrails [10,11], For this reason, much 
effort has been dedicated to developing physically repre- 
sentative ice cloud models (e.g., ice habits, particle size 
distributions, and mass-dimension relationships) and 
accurately computing the single-scattering properties of 
individual ice particles (i.e., the extinction efficiency, the 
single-scattering albedo, and the phase matrix). 

However, solving Maxwell's equations for light scatter- 
ing by ice particles from first principles is much more 
complex than light scattering by water droplets. Ice 
particles are nonspherical and morphologically complex, 
while having a broad range of particle sizes ranging from a 
few microns to centimeters. In a historical context, the 
majority of simulations [12-19] are based on geometric- 
optics principles that are asymptotically correct and 
approximately valid at large size parameters, i.e., when 
the particle size is sufficiently large relative to the wave- 
length. The ratio of the particle circumference to the 
wavelength is known as the size parameter. A number of 
physical-geometric optics methods have also been pro- 
posed to improve the accuracy in geometric-optics 
approximations [20-26], 

Deschamps [27] addressed the importance of the “ray” 
concept of light in engineering, stating that “if the nature 
of light and Maxwell's equations had been known, earlier 
optical instruments would not have been invented so 
readily"! Similarly, the application of ray techniques in 
solving light scattering by individual ice particles in the 
atmospheric radiation discipline plays a critical role 
because it provides a first-order approximation for higher 
frequency scattering where rigorous solutions are unavail- 
able. In other words, the knowledge of ice cloud radiation 
would have been severely limited in a modeling perspec- 
tive without employing the geometric-optics method. 
A more specific example is that for an ice particle with 
a size parameter larger than 200, the ability to obtain 
a rigorous set of optical properties is almost beyond our 
current first-principle modeling capabilities. For this rea- 
son, the geometric-optics methods continue to be popular, 


but there is a growing awareness of the need to provide an 
accuracy assessment due to the inherent approximations 
[e.g., 28, 29], 

With the development of computational electrody- 
namics, numerical techniques to solve Maxwell’s equa- 
tions are now available to obtain numerically exact 
solutions. The power of these numerical techniques began 
to have an impact within the past two decades. A few 
frequently used computational techniques are the finite- 
difference time-domain (FDTD) [30-33], the pseudo- 
spectral time-domain (PSTD) [34-36], the discrete- 
dipole-approximation (DDA) [37-40], and the T-matrix 
method [41-45]. Based on the numerical techniques, a 
much better knowledge has been obtained of the optical 
properties of particles at small-to-moderate size para- 
meters. The applicability of the geometric-optics in differ- 
ent size parameter domains is investigated whenever 
possible. For example, Mishchenko and Macke [46] inves- 
tigated T-matrix computations of light scattering by cir- 
cular cylinders and found that ice halos, an optical 
phenomenon predicted from geometric optics, would not 
be observed when the particle size parameter is less than 
approximately 100. Unlike the geometric-optics methods, 
the performance of numerical methods strongly depends 
on computational resources. As computer power increases, 
some new domains of rigorous solutions are being con- 
quered. Meanwhile, because of its simplicity, the use of the 
geometric-optics method is preferred for a complete range 
of ice particle sizes in the wavelength spectrum from UV to 
the near-infrared. 

To our knowledge, the induced uncertainties due to the 
nature of the approximations inherent in the ray-tracing 
technique have not yet been assessed in remote sensing 
and climate studies involving ice clouds, most likely 
because the rigorous solution domain is still extremely 
limited. The Invariant Imbedding T-matrix method (II-TM) 
[47-50] is applicable to a broad range of size parameters 
from the Rayleigh region up to the geometric optics 
domain where ice halos are observed. In this study, we 
use the II-TM to assess the accuracy of geometric-optics 
approximation and the resulting uncertainties in remote 
sensing and radiative transfer simulations. 

The theoretical components of the employed geometric- 
optics method are delineated before the assessment because 
a number of modifications in geometric-optics methods exist 
in the literature. A brief review of available geometric-optics 
methods is provided in Bi and Yang [51], The geometric 
optics methods most frequently used are the conventional 
geometric-optics method (CGOM) and the improved 
geometric-optics method (IGOM), although a few more 
rigorous, but relatively less computationally efficient, meth- 
ods have also been developed [20-26], One unique differ- 
ence between the CGOM and the IGOM is that the latter 
considers the spreading of scattered beams (a physical-optics 
effect) when propagating from the near-field to the far-field 
region. As the particle size parameter increases, the IGOM 
simulated phase matrix transitions to that from the 
CGOM when the spreading effect is negligible. With the 
ray-spreading effect incorporated, the IGOM is applicable to 
relatively small size parameters (20-100) where the ice halos 
disappear, but the geometric-optics principles are still more 
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or less valid in estimating the internal electromagnetic field. 
The CGOM and the 1GOM also differ in the computation of 
the extinction efficiency Q ex[ , defined by the ratio of particle 
extinction cross section to particle projected area. The CGOM 
empirically assumes that the extinction efficiency is constant 
at a value of 2 in conjunction with the blocking effect and the 
diffraction effect whereas the 1GOM computes the extinction 
efficiency based on a physical-optics approximation that 
takes into account the interference between the diffraction 
and the forward scattered rays. 

In this study, we begin with the assessment of the 
accuracy of the single-scattering properties derived from 
the CGOM and the IGOM in the T-matrix applicable 
domains. Next, we use a combination of the 1I-TM and 
the IGOM for computing the single-scattering properties in 
a complete range of particle sizes ranging from 2 to 
10,000 um in maximum diameter. Based on the finding 
that the IGOM and the T-matrix results are very similar at 
the transition regions from wave optics to geometric 
optics, we consider the Il-TM+IGOM results as a bench- 
mark and assess the uncertainties of applying the CGOM to 
an entire range of particle size parameters. From an 
application perspective, we then assess the uncertainty 
of using the CGOM in retrieving the ice particle effective 
diameter and optical thickness from Moderate Resolution 
Imaging Spectroradiometer (MOD1S) measurements and 
the broadband radiative forcing calculations using the 
general circulation model (GCM) version of the Rapid 
Radiative Transfer Model (RRTMG), which is the model 
employed in the Community Atmosphere Model (CAM, 
version 5.1) developed at National Center for Atmospheric 
Research. For simplicity, the ice particles are assumed to be 
smooth-faceted hexagonal columns. Although atmo- 
spheric ice habits are diverse, the accuracy of the CGOM, 
IGOM, and II-TM is not assumed to change for different 
habits. 

The remainder of this paper is organized as follows. 
In Section 2, we describe the numerical methods used 
for computing the optical properties of ice particles. The 
numerical results for the optical properties of hexagonal ice 
particles are shown in Section 3, including the comparison 
of the single-scattering properties and bulk-scattering prop- 
erties computed from the CGOM and the II-TM + IGOM. 
Section 4 provides sensitivity studies to assess the uncer- 
tainties of applying the CGOM in the applications of 
remote-sensing retrievals. In Section 5, the uncertainties 
in radiative forcing calculations due to applying the CGOM 
are quantified through a single-column radiative transfer 
model and a GCM. The conclusions of this study are given in 
Section 6. 


2. Computational methods 

In this section, we briefly outline the principles of the 
computational methods of the CGOM, the IGOM, and the 
II-TM. The readers are referred to [51,52] as primary 
references for the applications of geometric optics. 
A common theoretical basis of the CGOM is a combination 
of Fraunhofer diffraction by the particle shadow (Fig. la), 




Fig. 1 . Conceptual figure of geometric-optics principles. Pdiffaaion indi- 
cates the phase function associated with Fraunhofer diffraction. P ray C gom 
and Pray, igom are phase functions associated with the scattered rays 
computed from the CGOM and the IGOM, respectively. In the CGOM, 
the phase function P ray is computed from the division of energy A E with 
respect to the differential solid angle A £2. 


and the angular distribution of scattered intensity 
obtained from the ray-tracing technique shown in 
Fig. lb. Specifically, the scattering phase function from 
the ray-tracing process is obtained from the division of the 
energy received by each detector and the solid angle. The 
extinction efficiency is assumed to be 2 and arises sepa- 
rately from the diffraction and the ray-tracing part. The 
absorption efficiency is computed based on the energy lost 
in the ray-tracing process. The assumption of a constant 
value of Q ext in the CGOM may affect the accuracy of the 
single-scattering albedo, 1 - Q. a bs/Q.ext< where Q sca and Q abs 
are the scattering and absorption efficiencies, respectively. 
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The largest advantage of the CGOM is in its simplicity and 
flexibility to adapt to particle geometries; however, it 
suffers from several drawbacks. The solid angle is zero in 
both the forward and backscattering directions, thus 
causing a singularity problem. In most CGOM algorithms, 
values of the phase function in the direct forward and 
backward directions are obtained through extrapolation of 
the near-forward and near-backward directions. Moreover, 
if the particle is in a fixed orientation, the phase function 
is not a continuous curve but several discrete points. 
Although the random-orientation condition produces a 
relatively smooth phase function, an isolated peak is 
usually observed in the direction of forward scattering, 
commonly called the delta-transmission [15] energy, that 
is fundamentally related to parallel facets in a hexagonal 
ice habit. In the case of distorted or roughened particles, 
the delta-transmission term is not observed [53], 
Mishchenko and Macke [54] specifically explained that 
the delta-transmission term exists because of the missing 
physical-optics effect. When the particle is non-absorptive, 
the angular distribution of scattered energy is evidently 
size-independent, but this is contradictory to the wave 


Table 1 

Theoretical components involved in the CGOM and the IGOM. 



CGOM 

IGOM 

Extinction efficiency 

2 

Asymptotic value is 2 

Delta transmission 

Removed 

Removed 

Ray spreading 

No 

Yes 

Inhomogeneous wave 

Yes 

Yes 


theory of light. Furthermore, the extinction efficiency 
is an asymptotic value and does not capture the size- 
dependence for small and moderate size parameters. 

To overcome the shortcomings of the CGOM, 
geometric-optics principles are employed to compute 
either the internal field or the particle surface field, and 
the far field is obtained through exact electromagnetic 
volume integral equations or surface integral equations 
[20-26]. Within the new framework, the physical-optics 
effects are taken into account. For example, the extinction 
efficiency considers the interference between the diffrac- 
tion and the waves associated with the forward scattering 
rays. The diffraction effect associated with outgoing scat- 
tered beams (Fig. lc) is taken into account so that the 
scattering phase function from the ray-tracing calculation 
is size-dependent. For example, halos will only appear 
when the particle size parameter is sufficiently large (more 
than 100), where the ray-spreading effect is small and the 
rays are strongly localized. Yang and Liou [21] proposed a 
simplified algorithm to account for the physical-optics 
effect in the phase matrix computation that is implemen- 
ted in IGOM. 

When ice particles are absorptive, most localized waves 
inside the particle are inhomogeneous. In this study, we 
have considered the inhomogeneous nature of localized 
waves in the ray-tracing process in both the CGOM and 
the IGOM, although for simplicity the effect is usually 
neglected in many studies. Two obvious reasons to justify 
the consideration of inhomogeneous waves are that: (1) 
the energy is not conserved at the first-order of refraction 
if the inhomogeneity of localized waves is omitted; (2) the 
asymptotic value of the extinction efficiency derived from 








Fig. 2. Comparison of phase functions computed from the 11-TM, the IGOM and the CGOM. The wavelength and the hexagonal column height are indicated 
in each sub-figure. The size parameters for each sub-figure from (a)-(f) are approximately 126, 110, 118, 129, 121 and 124. 
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the IGOM differs from 2 by a non-negligible amount if not 
considering the electric field component along the propa- 
gation direction. 

In consideration of a number of geometric-optics algo- 
rithms, we specifically delineate the theoretical compo- 
nents involved in the CGOM and the IGOM in Table 1. Note 
that the delta-transmission peak has been removed in the 
CGOM by considering the physical-optics effect. 

For a particle having a small size parameter, where the 
ray-concept is inappropriate for calculating the electro- 
magnetic field inside the particle or particle surface, the 
IGOM gives inaccurate results; Maxwell's equations need 
to be solved. The transition (or T-) matrix method in 
conjunction with the extended-boundary condition (EBC) 
technique in the T-matrix computations pioneered by 
Waterman [41,42] is a powerful approach for computing 
the optical properties of nonspherical particles. In princi- 
ple, the EBC is applicable to an arbitrarily shaped non- 
spherical particle; however, the EBC is most successful in 
cases of axially symmetric particles. For example, for a 
hexagonal particle whose geometry is of 6-fold rotational 
symmetry, the EBC's applicable size parameter region is 
substantially narrowed in comparison to spheroid cases. 


Note that in the EBC method, the T-matrix is computed 
from surface integrals over the particle surface. Alterna- 
tively, the T-matrix can be derived from an electromag- 
netic volume integral equation and computed from an 
invariant imbedding procedure. One unique advantage of 
the invariant imbedding approach is that the algorithm is 
stable and flexible in adapting to a given particle geometry. 
We have extended the II-TM applicable region to large size 
parameters for nonspherical particles [47-49], Here, we 
will not discuss the theoretical basis of the II-TM as the 
details can be found in [47-50], 


3. Optical properties of ice particles 

The optical properties of hexagonal ice particles are 
computed in a spectral region from 0.2 to 100 pm. The 
aspect ratio of hexagonal ice particles is the same as the 
geometry defined in [55] (their Table 1). Six wavelengths 
in the visible and infrared are selected for remote-sensing 
studies: 0.65, 0.86, 2.13, 8.5, 11, and 12 pm. The ice 
refractive indices compiled by Warren and Brandt [56] 
are used in the simulations. 







Fig. 3. Comparison of the extinction efficiency, the single-scattering albedo, and the asymmetry factor at wavelengths (a) 0.65 pm (left column), (b) 
2.13 pm (middle column), and (c) 12 pm (right column) computed from the II-TM + IGOM and the CGOM. 
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Fig. 2 shows a comparison of phase functions computed 
from the CGOM, the II-TM and the IGOM at six representa- 
tive wavelengths. The particle sizes for each phase function 
are indicated in the figure. The phase functions computed 
from the II-TM and IGOM are very similar. At wavelengths 
of 0.65 and 0.86 pm (Figs. 2(a) and (b)), the ice particles are 
non-absorptive and ice halos at 22° and 46 are beginning 
to emerge for the selected size in the II-TM and IGOM 
calculations, but they are much too prominent in the CGOM 
results. Flowever, a peak that is predicted by the IGOM and 
the CGOM at approximately 150° is not observed in the II- 
TM simulated phase function. The size parameter is most 
likely insufficiently large, so that the ray-concept in the 
CGOM and the IGOM fails for higher-order rays. Addition- 
ally, a small peak of phase function at the wavelength of 
8.5 |jm predicted from the CGOM is not observed from the 
II-TM simulation. Because the II-TM is computationally 
expensive and the II-TM and the IGOM results are very 
close as shown in Fig. 2, we only employ the IGOM in the 
following simulations when the particle sizes are larger 
than the numbers indicated in the figure. The IGOM results 
of the extinction efficiency and the absorption efficiency are 
adjusted with respect to the II-TM for a smooth transition 
through a semi-empirical estimation of the edge effect 
contribution [46], 

Fig. 3 shows a comparison of the extinction efficiency 
(Qext). the single-scattering albedo ( a > ), and the asymmetry 
factor (g) of hexagonal columns at three wavelengths (0.65, 
2.13, and 12 pm), computed from the CGOM and II- 
TM + IGOM, and for particle sizes ranging from 2 to 
10,000 pm. The differences between the results from the 
two techniques are evident, particularly for small particle 
sizes. Note that the extinction efficiencies and the asym- 
metry factors at the wavelengths of 0.65 and 2.13 pm 
oscillate, but the CGOM results do not because the inter- 
ference effect is not taken into account. While the single 
scattering albedo is well approximated by the CGOM at 
2.13 pm, large errors for small particles occur at 12 pm. Note 
that the size parameter for a given particle size at 12 pm is 
approximately 5.6 times smaller than that at 2.13 pm and 
the particle is more absorptive at 12 pm than at 2.13 pm. As 
is evident from the comparison, the T-matrix and IGOM 
solutions cover a size parameter range where the conven- 
tional geometric-optics principles break down. 

Based on the independent scattering assumption, ice 
cloud bulk-scattering properties are obtained by integrat- 
ing the single-scattering properties of ice particles with 
the particle-size distributions (PSD). The Gamma size 
distribution is adopted in this study, and the PSD can be 
represented by the following function: 

n(L) = N 0 (L/v) a ~ A exp( — L/p), (1) 

where N 0 is a normalization factor and L is the maximum 
dimension (length of the column). The PSD is determined 
by a and u, the shape and scale parameters. To determine 
the two parameters, we use the values from ice cloud 
in situ measurements [57] and from retrievals [58] for 
contrails. Fig. 4 shows the relationship between a and v. 
The cirrus PSD data compiled by Baum et al. [57] are 
illustrated by blue dots in the figure, and the retrieved data 
from Iwabuchi et al. [58] are shown as green circles. The 



Baum et al. 201 1 [57] 
o Iwabuchi et al. 201 2 [58] 
— Eq. (2) 


10 '' 10 ° 10 ' 10 2 10 3 
Scale parameter v (pm) 

Fig. 4. Size distributions used in the computation of ice cloud bulk- 
scattering properties. The red dashed curves in the figure indicate that 
most data lay in the range either twofold smaller or larger than the fitting 
curve. (For interpretation of the references to color in this figure legend, 
the reader is referred to the web version of this article.) 


data show a large dispersion. Different from the linear fit 
given by Iwabuchi et al. [58], this study uses a hyperbolic 
tangent function to fit the observations, and the solid red 
curve in the figure follows: 

log(a+2) = 0.44-0.32 x tanh [2 x log(v)-4], (2) 


A total of 18 pairs of a and v values corresponding to the 
effective particle diameter D eJ y ranging from 10 to 180 pm 
at an interval of 10 pm are used in this study. D e ff is defined 
by 


D 


3<V) 

eff 2(A) ’ 


(3) 


where (17) and (A) are the averaged volume and projected 
area of the particles associated with a given size distribu- 
tion. The dashed curves in the figure indicate that most 
data lay in the range either twofold smaller or larger than 
the fitting curve. 

Fig. 5 compares the bulk-scattering phase functions 
computed from the CGOM and the II-TM + IGOM at 0.65, 
2.13, and 12.0 pm. The differences are evident when the 
D e ff is small, but become negligible when D e ff becomes 
sufficiently large. Indicated in the figure are the bulk 
extinction efficiency, the bulk single-scattering albedo, 
and the bulk asymmetry factor. The edge-effect con- 
tribution to the extinction efficiency is observable. In 
the case of a large D e ff, small particles contribute less to 
the total cross section because of their smaller pro- 
jected area and relatively low weight of probability in 
the particle size distribution. The difference between 
the extinction/scattering cross-sections computed by 
the II-TM + IGOM and the CGOM decrease as the D ef f 
increases. 

From the comparison of the single-scattering proper- 
ties computed from the CGOM and the “benchmark”, we 
identify errors in optical properties of small ice crystals 
calculated using CGOM. In the following sections, we 
quantify the radiative uncertainties if the CGOM is applied 
to all particle size parameters. 
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a 




b 




c 






Scattering Angle ( 0 ) 


Scattering Angle ( 0 ) Scattering Angle ( 0 ) 


Fig. 5. Comparison of bulk ice cloud scattering phase functions (P n ) for different effective diameters (D e ff) at the wavelengths of (a) 0.65 pm (left column), 
(b) 2.13 pm (middle column), and (c) 12.0 pm (right column). Indicated in the figure are the bulk extinction efficiency (Q exf ), the bulk single-scattering 
albedo (a), and the asymmetry factor (g) of the bulk-scattering phase function. 


4. Remote sensing implications 

The previous section compared both the single- and 
bulk-scattering properties of hexagonal columns com- 
puted from the II-TM + IGOM and the CGOM. The effect 
of the bulk-scattering property differences on the remote 
sensing implementations is assessed in this section. Two 
independent retrieval methods are used that employ the 
hexagonal column model and the scattering properties 
obtained from the II-TM + IGOM and the CGOM. One 
retrieval method is based on two solar reflectance bands 


and the other on three infrared (IR) bands. The accuracy of 
the CGOM is evaluated by comparing the retrieved ice 
cloud properties (i.e„ cloud optical thickness x and D eJ y) 
with those obtained from the II-TM+IGOM simulations for 
these two methods. 

Both retrieval methods use measurements from the 
Moderate Resolution Imaging Spectroradiometer (MODIS) 
on the National Aeronautics and Space Administration 
(NASA) Terra and Aqua platforms. MODIS has 36 spectral 
bands ranging from 0.41 to 14.2 pm. In the current study, 
we employ the Aqua/MODIS level-1 B (LIB) collection-5 
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(“MYD021KM”) product that provides top of the atmo- 
sphere (TOA) radiances/reflectivities for the 0.86- and 
2.13-pm bands and the brightness temperature for three 
IR bands (8.5-, 11- and 12-pm), which are used in the 
retrieval procedures, and a product ("MOD03”) that con- 
tains the geolocation and viewing geometry information 
for the satellite observations. The MODIS-retrieved cloud 
phase, optical thickness, and top height from the Aqua/ 
MODIS level-2 cloud product (MYD06) are also used. The 
MOD1S pixels that are over the ocean and identified as ice 
clouds in the MODIS level-2 products are used for our 
retrievals. For the IR retrieval, which is inherently less 
sensitive to optically thick clouds, the cases with MODIS 
optical thicknesses larger than 6 are neglected. The atmo- 
spheric profiles used in the RT calculations are provided by 
the Modern Era Retrospective-Analysis for Research and 
Applications (MERRA) model, specifically the 3-hourly 
profiles (“Int3_3d_ams_CP”) of temperature, water vapor 
density, and ozone density defined for 42 pressure levels 
at a grid resolution of 1.25° x 1.25°. 

The first retrieval algorithm is based on two solar 
reflectance bands: specifically, a weakly absorbing, visible 
or near-infrared window band (VIS/NIR) (e.g., 0.64 or 
0.86 pm) that is sensitive mainly to r, and an ice absorbing 
shortwave infrared (SW1R) band (e.g., 1.6 or 2.13 pm) that 
is sensitive to both r and D e ff [59], This study uses the 
MODIS observations at 0.86 and 2.13 pm bands, and the 
fast radiative transfer model (RTM) developed by Wang 
et al. [60] to simulate the reflectance of the modeled ice 
clouds. The fast RTM uses pre-computed bidirectional 
reflectance/transmittance distribution functions (BRDF/ 
BTDF) for single layer clouds, and lookup tables (LUTs) 
for different r and D eJ y are calculated using the discrete 
ordinates radiative transfer program (DISORT). Fig. 6 
shows the comparisons of the LUTs built on the optical 
properties calculated by the II-TM+IGOM and the CGOM, 
respectively. The solar zenith angle is assumed to be 20°. 
Two viewing zenith angles, 37° and 70°, are selected for 
radiative transfer simulations to build LUTs as shown in 
Fig. 6a and b, respectively. Substantial differences exist in 
the D e ff isolines, because the differences between the 


optical properties given by the II-TM + IGOM and the 
CGOM become larger as D ef f decreases. Furthermore, as 
shown in Fig. 6, the relative differences between the LUTs 
are sensitive to the solar-satellite geometry, because the 
relative magnitude between the CGOM and the II- 
TM+IGOM phase functions change with respect to the 
scattering angles (see Fig. 5). For the case in Fig. 6a, the D e ff 
isolines computed from the CGOM optical properties are 
lower than the II-TM+IGOM counterparts, which implies 
that the retrieved D ej j based on the CGOM in all probability 
be smaller than those based on the II-TM+IGOM. How- 
ever, for the case in Fig. 6b, the CGOM optical properties 
result in higher reflectivities at the same D e ff, suggesting 
that the retrieved D e jj is mostly likely to be overestimated. 
The two LUTs of the CGOM and the II-TM+IGOM in Fig. 6a 
and b converge as D e ff reaches 100 pm and 30 pm, respec- 
tively. In terms of the r isolines shown in Fig. 6, the results 
for the II-TM + IGOM and CGOM are similar with notice- 
able differences for small D e ff. 

To retrieve the cloud optical thickness and effective 
particle diameter, the following cost function is defined 

2 (K 0.86 —Rs, 0.86+ , (Ro, 2.13 ~Rs, 2.13+ ... 

X = { Ro, 0.86 J + { Ro, 2.13 J - (4) 

where R 0 x and R s , a are the MODIS bidirectional reflectiv- 
ities and the fast RTM simulations at wavelength A, respec- 
tively. The r and D e ff values are obtained from the simulated 
reflectivities that minimize Eq. (4). Fig. 7 illustrates the 
cloud properties retrieved from the MODIS measurements. 
A MODIS granule at 07:45 UTC on August 2, 2010 is used in 
the retrieval. Fig. 7a and b shows the solar zenith and 
viewing zenith angles, respectively. Fig. 7c and d provides 
the retrieved cloud r and D e j based on the LUT developed 
from the scattering properties given by the II-TM+IGOM. 
As illustrated from the retrieved results, the clouds in this 
granule have a wide range of r and D,,jf values. 

The relative differences between the II-TM+IGOM and 
CGOM cases are shown in Fig. 7e and f, and the relative 
differences are given by 

RE t = *CCOM-*U-TM + ICOM x 1()0% (5) 

T II-TM + IGOM 



Fig. 6. Comparisons of the LUTs using 0.86- and 2.13-pm reflectance calculated using scattering properties from the II-TM + IGOM and the CGOM. 
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Fig. 7. The solar zenith (a) and view zenith (b) for the MODIS observations, the retrieved cloud optical thickness (c), and the retrieved effective diameter 
(d) based on the II-IM + IGOM, and the relative differences of the retrieved optical thickness (e) and effective diameter (f) based on the II-TM + IGOM and 
the CGOM databases. 


and 


RE n 


E>eff,CGOM - %/,» - TM + /COM x 10Q% 
DeffJI-TM + lGOM 


( 6 ) 


where th_ tm+ i CO m and D e ffj I _ TM+IC0 M, and t C com and 
Deff.ccoM are the values retrieved with the scattering 
properties given by the II-TM + IGOM and the CGOM, 
respectively. As shown in Fig. 7e and f, the relative 
differences caused by the CGOM are as large as 15% and 
30% for z and D r ,jj, respectively. To be more specific, the 
relative differences of t, generally under 10%, show weak 
dependence on the satellite viewing geometry, and the 
largest differences occur at those pixels with small t. 
However, it should be noted that the relative differences 
on the retrieved D eJ y shown in Fig. 7f are much larger, and 
strongly depend on the viewing zenith angle. The 


correlation between Fig. 7b and f can be explained by 
the LUTs shown in Fig. 6. The CGOM-based D e // isolines in 
the LUTs can be lower or higher than the II-TM+IGOM 
counterparts because the reflectivity for a specific solar- 
satellite configuration is highly related to the magnitude of 
the phase functions at the corresponding scattering angles 
where the CGOM phase function can be larger or smaller 
than the II-TM+IGOM phase function. Consequently, the 
differences of retrieved effective diameters from the two 
LUTs can be smaller or larger with a similar pattern as the 
viewing zenith angle. 

Fig. 8 shows the mean biases (solid line) and the 
standard deviations (shaded areas) of the retrieved r and 
D e .ff values for the CGOM-derived scattering properties. It 
can be seen that the retrieved t is quite accurate for t < 30 
with mean biases of less than 1% and standard deviations 
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Cloud Optical Thickness 



Fig. 8. The mean and standard deviations for the relative difference in the retrieved optical thickness (a) and effective diameter (b). 


less than 5%, but the CGOM-based retrievals underesti- 
mate z by approximately 5% and have standard deviations 
that can exceed 5% as z increases. Note that, the reflectance 
at 0.86 pm becomes less sensitive to r for thicker clouds. 
For example, the reflectivity increases by an amount < 0.1 
as the t increases from 30 to 100 for all D e ff. Low sensitivity 
indicates the existence of large retrieval differences for 
thick clouds. Meanwhile, the relative differences for the 
retrieved D e jj are more significant with mean biases and 
standard deviations of approximately up to 20% and 40%, 
respectively, for the particles with D eJ y< 50 pm. However, 
the errors become much smaller as D e ff increases. 

Fig. 9 shows the relative differences in CGOM-retrieved 
D e ff within different ranges of z. As evident from the figure, 
the errors of D e jy are very sensitive to z. For example, the 
patterns in Fig. 9a and b are quite different and very large 
bias are obtained for optically thin clouds with small 
particle sizes. Note that the pixels for D ej y> ~85 are not 
available when r< 1. In comparison of Fig. 9c and d, the 
errors of D e // are relatively larger for z < 2 than those for 
z > 2. As shown in Fig. 9e and f, the errors in D e ff are most 
significant for cases where z < 4, and much smaller as 
z > 4. Because Fig. 9e is similar to Fig. 8b, the large errors 
shown in Fig. 8b primarily occur when r < 4. Different 
from the above findings, the retrieved z is found not to be 
very sensitive to D eJ y. For example, the pattern of the 
relative differences of z are very similar for the pixels with 
D e ff< 50 pm as for those with D e g> 50 pm (figure not 
shown). 

The second retrieval method is based on the three MODIS 
1R bands at 8.5, 11, and 12 pm. The fast high-spectral- 
resolution RTM (HRTM) developed by Wang et al. [61,62] is 
used to simulate the radiances and resulting brightness 
temperatures for the three bands. Gaseous absorption is 
taken into account by a pre-computed transmittance data- 
base using a rigorous line-by-line radiative transfer model 
(LBLRTM), and the scattering properties given by the II- 
TM+IGOM or the CGOM are used to calculate the angular- 


dependent ice cloud reflectance and transmittance, the 
effective emissivity, and the effective temperature functions. 
The spectral response functions of the three MODIS bands 
are considered in the HRTM. The t-D^ pair that leads to the 
least root mean square (RMS) value is chosen as the most 
appropriate inference of the ice cloud properties, and the 
RMS is defined as 


RMS = 


(BT 0 gs -BT S 8 . 5 ) + (BT 0 11 — BT S n) +(BT 0 l2 - BT Si i2 Y 


1/2 


( 7 ) 


where BT 0i 2 and BT S> 2 are the brightness temperatures from 
MODIS observations and the HRTM simulations at wave- 
length X. Fig. 10 is the same as Fig. 8 but for z and D e j from 
the IR retrieval. The retrievals carried out using the scattering 
properties from the CGOM systemically overestimate the r by 
approximately 12% with standard deviations over 3%. For D ej f, 
the mean biases caused by the inaccuracy of the CGOM range 
between about - 10% and 10% with standard deviations up 
to 20%. 


5. Uncertainties in radiative forcing calculations 

The differences in the ice cloud optical properties 
between the II-TM+1GOM and CGOM cases also influence 
the parameterization schemes of ice cloud bulk optical 
properties used in the single-column radiative transfer 
models (RTMs), as well as in general circulation models 
(GCMs), thereby affecting the simulations of ice cloud 
broadband radiances. Results are shown for bands in the 
shortwave (SW) spectrum (X < 5 pm) and in the longwave 
(LW) spectrum (5 < X < 100 pm). The cloud radiative effect 
(CRE) is defined as the difference in the net TOA/surface 
flux (SW, LW or total) between the cloudy sky and the 
clear sky. 
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Fig. 9. The mean and standard deviations for the relative difference in the retrieved effective diameter within different ranges of the optical thickness. 
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Fig. 10. The mean and standard deviations for the relative difference in the retrieved optical thickness (a) and effective diameter (b) for the 1R retrieval. 


The optical properties of ice clouds used in numerical 
model simulations include the mass extinction coefficient 
( k ext ), the mass absorption coefficient k abs , the bulk single- 
scattering albedo m, and the mean asymmetry factor g, which 


are defined by 

It ft™ Q_ ex t(L,A)A(L)S(A)n(L)dLdA 
ft ft V(L)S(A)n(L)dLdA 


( 8 ) 
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Fig. 11. Comparison of bulk ice cloud optical properties derived from 


k irr JiT ^(L, A)A(L)S(A)n(L)dLdA 

itr /i”T V(L)S(A)n(L)dLdA 

( 9 ) 

co= \ — k abs / k ex t , 

( 10 ) 

itr fi;™g<.CA)Q. SC a(L,A)A(L)S(A)n(L)dLdA 
iir Jr™ Q-sca(L,A)A(L)S(A)n(L)dLdA 

( 11 ) 


where L is the height of the hexagonal column in a particle 
size range [Lmin, f-max]. V and A are the volume and the 
average projected area, /l^n and A m ax specify the wavelength 
limits of a given SW/LW band, and S(2) is the solar spectrum 
for SW bands and the Planck function (the cloud temperature 
is assumed to be 233 K) for LW bands. All are parameterized 
as functions of ice cloud D e g. Fig. 11 shows the ice cloud bulk- 
scattering properties for selected spectral bands of the RRTMG 
radiative transfer model. For the CGOM case, the mass 
extinction coefficient is invariant with the SW spectral bands, 
and will cause a difference in the mass extinction coefficient 
of up to 0.05 m 2 g~ ] (14%) for the 2.15-2.5 pm band for small 
D e ff. Similarly, the mass absorption coefficient for a given LW 
band (e.g., the 14.3-15.9 pm band) can differ by about 
0.04 m 2 g _1 (24%). The largest differences in the asymmetry 


b 
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d 



Effective diameter (pm) 


l-TM+IGOM and CGOM at selected spectral bands of the RRTMG RTM. 

factor occur for small D e j in the longer wavelength bands. 
Relatively small differences are found for the single-scattering 
albedo (SSA), where the I1-TM+1GOM case has slightly lower 
SSA than the CGOM case for the near IR band. 

Although differences in optical properties are evident, 
it is useful to understand how these differences affect the 
simulated ice cloud CRE, especially for single bands in the 
SW and LW spectrum. Here, the RRTMG SW/LW model 
[64] is implemented with a standard mid-latitude summer 
atmospheric profile, assuming a solar zenith angle of 60°, a 
single-layer ice cloud located at 12 km, a surface albedo of 
0, and an emissivity of 1 for all wavelengths. The left panel 
in Fig. 12 shows the differences of TOA SW, LW, and total 
(SW+LW) ice CRE obtained from CGOM and II-TM + IGOM 
databases as functions of D e //and % due to the difference in 
optical properties. The difference is obtained by subtract- 
ing the CGOM-based CRE from the II-TM+lGOM-based 
CRE. The SW radiative effect is negative, meaning 
that II-TM + IGOM case results in a stronger ice cloud SW 
radiative effect than the CGOM case. The largest differ- 
ences can be up to 7 W m~ 2 for r ranging from 5 to 10 and 
D e ff =10 pm (i.e„ small values). The LW radiative effect 
difference is positive, and twice as large as that of the SW 
at small D eJ ^ when t is near 1. The resulting total CRE 
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Fig. 12. Left panel: RRTMG RTM simulated TOA cloud radiative effect difference between the Il-TM+IGOM and CGOM cases: solar band CRE (a), IR band 
CRE (b), and total (solar+lR) CRE (c). Right panel: the percentage of the relative difference in comparison with the II-TM+IGOM reference. 


differences are dominated by the LW feature, resulting in 
positive differences for small optical thickness values 
(r<5) but negative differences for large values of optical 
thickness (r > 5). The right panel shows relevant difference 
percentages. Although the differences of TOA SW tend to 
be small in the cases of low t and D e ff values, the relevant 
percentage differences may be significant. A similar fea- 
ture is observed in the LW case. 

To identify the contribution of CRE from various spec- 
tral bands of RRTMG RTM, the band-by-band CRE is 
calculated. Fig. 13 shows four SW/LW spectral bands with 
the largest contributions to the CRE difference between 
the II-TM + IGOM and the CGOM case. The corresponding 
differences in the optical properties are clearly shown in 
the CRE at the spectral bands. 


Following Yi et al. [65], the bulk ice cloud optical 
properties derived from the II-TM+IGOM and the CGOM 
cases are implemented as parameterizations into the NCAR 
Community Atmospheric Model (CAM, version 5.1 ). A ten- 
year climatology is derived from the model runs with the 
two parameterizations, and the global annual averaged 
CRE are calculated and shown in Fig. 14. The SWCRE 
(Fig. 14a) is negative; whereas, the SWCRE difference 
(Fig. 14d) is generally positive, indicating that the 
II-TM+IGOM case results in a larger SWCRE. At the same 
time, as shown in Fig. 14b and e, the II-TM+IGOM case 
also has a stronger LWCRE. While the SW/LW CRE differ- 
ence can be up to 8 W m -2 regionally, for example, in the 
tropical area, the global averaged SW and LW CRE differ- 
ences are 0.44 Wm -2 and — 1.35 Wm~ 2 , respectively. 
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Fig. 13. Band-by-band CRE differences between the II-TM+IGOM and the CG0M cases simulated by the RRTMG RTM at the top of the atmosphere, (a) TOA 
CRE ( Wm - 2 ), 0.78-1.21 [im, (b) TOA CRE (Wm~ 2 ), 0.63-0.78,1 m, (c) TOA CRE ( Wm ~ 2 ), 12.2-14.3,im and (d) TOA CRE ( Wm ~ 2 ), 10.2-12.2,im. 


The II-TM + IGOM-based net CRE (NTCRE) and the differ- 
ences between the CGOM and the II-TM + IGOM cases are 
shown in Fig. 14c and f, respectively. 

Based on the simulated results, the biases caused by the 
CGOM in the LW are significantly larger than in the SW. 
The physical reason is that for the longwave spectrum, the 
particle size parameters for a relatively large range of 
particle sizes are small so that the accuracy of the CGOM 
is rather poor. 

6. Summary and conclusions 

The numerical techniques based on geometric-optics 
principles and its modifications provide approximate solu- 
tions for light scattering by ice crystals. In reality, the 
geometric-optics approximation is widely used but gen- 
erally without an awareness of its uncertainties. Here, 
following [51], the geometric-optics methods are classified 
into the conventional geometric optics method (CGOM) 
and the improved geometric optics method (IGOM). 
A combination of the IGOM and the invariant imbedding 
T-matrix method (II-TM) is employed to derive the “bench- 
mark” single-scattering properties of hexagonal particles 
with diameters from 2 to 10,000 pm over a spectral range 
from 0.2 to 100 pm. Based on the single-scattering proper- 
ties, we compute the ice cloud bulk-scattering properties 
by using the particle size distributions obtained from 


in situ measurements. We begin with an assessment of 
the uncertainties of the CGOM in the computations of 
single- and bulk-scattering properties. Significant biases 
are identified in the optical properties of small individual 
ice particles calculated by the CGOM. Ice cloud bulk- 
scattering properties are applied in MODIS satellite retrie- 
val and ice cloud forcing simulations. In the MODIS 
satellite retrievals using the visible and shortwave infra- 
red, mean biases in optical thickness caused by applying 
the CGOM are found to be up to 5% with standard 
deviations of about 5%. Mean biases and standard devia- 
tions for the retrieved effective particle diameters can be 
over 20% and 30%, respectively, and generally decrease 
with particle size. However, mean biases in effective 
diameters for clouds with optical thicknesses above 4 are 
found to be below about 5%. For retrievals using IR bands, 
use of CGOM leads to systematic overestimation of the 
optical thickness by approximately 12% and biases in 
effective particle diameter ranging from - 10% to 10%. In 
the ice cloud radiative forcing calculations, biases are also 
found, in particular, from longwave radiative effect where 
the difference can be as high as 14 W m~ 2 , although global 
averaged SW and LW CRE differences are 0.44 W m~ 2 and 
-1.35 W m~ 2 , respectively. 

The sources of uncertainties in the retrieval and climate 
forcing calculation are primarily two-fold: the microphy- 
sical model and the accuracy of the computational 
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Fig. 14. Left panel: the Il-TM+IGOM-based ten year mean annual CRE (unit: W m-2): SWCRE (a), LWCRE (b), and NTCRE (c). Right panel: the differences 
between the I1-TM+1GOM and the CGOM cases (the CGOM case minus the ll-TM+IGOM case). 


techniques. In the present study, we only focus on the 
latter by assuming a simplified ice habit (i.e., smooth 
hexagonal column). In recent studies [63,66-68], the 
relative differences in retrieved optical thickness and 
effective particle diameter for different ice models can be 
comparable or larger than the maximum values found 
here. Similarly, in CRE calculations, differences in short- 
wave fluxes resulting from different microphysical models 
can be on the order of 10-20 W/m 2 [53,65] for ice particle 
roughness, and -70 to -30 [53] and -100 to -30 W/m 2 
[53,69] for non-unity particle aspect ratios. Actually, the 
uncertainties from 

the aforementioned two sources are not independent. 
For example, in the analysis of model uncertainties, the 
techniques used in the calculation of the single-scattering 
properties are assumed to be accurate. Further analysis of 
the uncertainties will be performed in future studies. 
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